浅谈C++浮点数精度的一些问题

有时,浮点数的精度是在做题中的一个重要的问题。举一个例子:

#include <stdio.h>
main() {
    printf("%.6lf\n", floor(sqrt(9999.9999) * sqrt(10000.0001)));
}

这份代码是一个简单的根号的乘法,我们可以手算得到 9999.9999×10000.0001=1000020.00012<10000\sqrt{9999.9999 \times 10000.0001} = \sqrt{10000^2 - 0.0001^2} < 10000,应当输出 9999.000000 。但是当你真正运行的时候,会发现它会输出 10000.000000 ,而不是 9999.000000 。这是为什么呢?这就要从浮点数的精度说起。

1. 浮点数精度和范围

在说浮点数精度之前,我们先说说一个相关的东西:科学计数法

用科学计数法表示的数一般是这样子的:符号a×bc\text{符号} a \times b^c ,其中符号为正或负数,aa 是有效数字,bb 是底数,cc 是指数。而且可以很容易知道,1a<101 \le a < 10 。例如在 3.14×1015-3.14 \times 10^{15} 中,- 是符号,3.14 是有效数字,10 是底数,15 是指数。

在计算机中,浮点数也用科学计数法表示。不过,因为用的是二进制存储,所以计算机中的浮点数一般这样表示:

符号a×2c\text{符号} a \times 2^c

例如,二进制下的 10101.10110101.101 (对应十进制下的 21.62521.625)就可以用科学计数法表示为(1.0101101)2×24(1.0101101)_2 \times 2^4

在C++中, float 是单精度,32位,double 是双精度,64位。并且这两个在内存中的分配结构也不一样:

符号位 指数 尾数
float 1bit 8bit 23bit
double 1bit 11bit 52bit

这里以 float 为例,符号位占 11 位,00 表示正数,11 表示负数。

指数占 88 位,表示 02550 \sim 255 ,为了能表示负指数,实际指数位存储在指数域中值减去一个偏移量(float127double1023),例如用 0000 0000 表示 0-127=1271000 0000 表示 128-127=11001 0001 表示 145-127=180001 0001 表示 17-127=-110 ,以此类推。(采用了移位存储)

尾数,也就是科学计数法中的底数,占 2323 位。显然底数的整数部分必定是 11 ,也就是这个小数必定是 1.×2c1. \cdots \times 2^c ,所以记的时候就把前面的 11 省略掉了,只记后面的小数部分。

可以发现,由于 223=83886082^{23} = 8388608 ,一共 77 位,这意味着最多能有 88 位有效数字(还要算上前面的 11),但绝对能保证的为 77 位,也就是说 float 的精度为 787 \sim 8 位有效数字;并且表示范围为 [2128,2128)[-2^{128}, 2^{128}) ,也就是 3.40×10383.40×1038-3.40 \times 10^{38} \sim 3.40 \times 10^{38}

同理,由于 252=45035996273704962^{52} = 4503599627370496 ,一共 1616 位,这意味着最多能有 1717 位有效数字(还要算上前面的 11),但绝对能保证的为 1616 位,也就是说 float 的精度为 161716 \sim 17 位有效数字;并且表示范围为 [21024,21024)[-2^{1024}, 2^{1024}) ,也就是 1.80×103081.80×10308-1.80 \times 10^{308} \sim 1.80 \times 10^{308}

举个例子:(12.625)10=(1100.101)2=(1.100101)2×(23)10(12.625)_{10} = (1100.101)_2 = (1.100101)_2 \times (2^3)_{10} ,所以这个浮点数在双精度下的存储是:

符号位 指数 尾数
00 1000 00101000 \text{ } 0010 1001 0100 000000001001 \text{ } 0100 \text{ } 0000 \cdots 0000100101100101 后面共 464600

2. 进制转换

整数的进制转换相信大家都会,那么这里只讲小数的进制转换。这里,小数的进制转换你会发现有三种情况:

(1) 有限小数变成有限小数 ,举个最简单的例子:(12.5)10=(30.2)4(12.5)_{10} = (30.2)_{4}

(2) 有限小数变成无限小数(或反过来) ,举个例子:(1.2)10=(1.00110011)2(1.2)_{10} = (1.00110011 \cdots)_2(103)10=(10.1)3(\frac{10}{3})_{10} = (10.1)_3

(3) 无限小数变成无限小数,举个例子:(103)10=(10.010101)2(\frac{10}{3})_{10} = (10.010101\cdots)_{2}

通过这些例子,可以发现,有些小数是无限的。

3. 精度失真

有些小数是无限的,但是精度是有限的,就会导致经过处理后的数据并不准确。这里我们举个例子方便理解:

(2.2)10=(10.00110011)2=(1.000110011)2×(21)10(2.2)_{10} = (10.00110011 \cdots)_2 = (1.000110011 \cdots)_2 \times (2^1)_{10} ,由于精度有限,转换后如下表所示:

符号位 指数 尾数
00 1000 00001000 \text{ } 0000 0001 1001 10011001 10100001 \text{ } 1001 \text{ } 1001 \cdots 1001 \text{ } 1010(共 1100010001111110011001 以及 1110101010,最后一位因为下一位为 11 所以进行了类似四舍五入的操作)

而上表中的数转换成十进制之后已经不是 2.22.2 了,而是 2.2000000000000002002.200000000000000200 。可以通过以下代码进行测试:

#include <stdio.h>
main() {
    printf("%.18lf\n", 2.2);
}

4. 例子的解释

我们来看一开始举的例子:

#include <stdio.h>
main() {
    printf("%.6lf\n", floor(sqrt(9999.9999) * sqrt(10000.0001)));
}

通过上面的解释,可以发现:

9999.99999999.99990000000079999.9999 \rightarrow 9999.9999000000007

10000.000110000.00009999999910000.0001 \rightarrow 10000.000099999999

开完算数平方根后(左边是失真的数开根后实际的数,右边是失真后的数):

99.99999950000000228500001142599.999999500000001299.999999500000002285000011425 \rightarrow 99.9999995000000012

100.000000499999995200000024100.0000004999999987100.000000499999995200000024 \rightarrow 100.0000004999999987

相乘后得到(左边是失真的数相乘后实际的数,右边是失真后的数):

9999.999999999999850010000.09999.9999999999998500 \rightarrow 10000.0

于是,原本到不了 1000010000 的数变到了 1000010000

5. 解决方法

其中一种解决的方法就是少用除法少开根号 。假如你要比较 ab\frac{a}{b}cd\frac{c}{d} 的大小,你就可以通过比较 a×da \times dc×bc \times b 的大小来减小失真(可以使得原来是有限小数的数(二进制下)还是有限小数,失真的概率小些)。假如你要比较 x2+y2\sqrt{x^2 + y^2}z (z>0)z \text{ } (z > 0) 的大小,你就可以直接比较 x2+y2x^2 + y^2z2z^2 的大小来减小失真。

还有一种解决方法就是手写浮点数,像这样:

...
struct BigDecimal {
    int before, after; //  分别表示小数点前的十进制数 和 小数点后的十进制数
    int pm = 0; //表示正负号, 这里 0 为正, 1 为负
    int length = 30, base = 1e8 + 0.5; // 表示小数点的长度 和 小数点后的进位制 (这里是 after 满 10^8 向前进位)
    BigDecimal(int pp = 0, int bb = 0, int aa = 0) {
        pm = pp; before = bb; after = aa;
    }
    bool operator < (const BigDecimal &A) const {
        if (pm != A.pm) return pm > A.pm;
        if (pm == 0) return before != A.before ? before < A.before : after < A.after;
        return before != A.before ? before > A.before : after > A.after;
    }
    bool operator > (const BigDecimal &A) const {
        return A < (*this);
    }
    bool operator == (const BigDecimal &A) const {
        return !((*this) < A || A < (*this));
    }
    ...
};
BigDecimal operator + (const BigDecimal A, const BigDecimal B) {
    if (A.pm == B.pm) return BigDecimal(A.pm, A.before + B.before + (A.after + B.after) / A.base, (A.after + B.after) % A.base);
    else ...
}
...

这里的原理就是本来小数点后的数会失真,将其变为整数后就不会失真了。

或者有的时候这样写可以更加简单(这里假设保留 44 位小数):

...
typedef long long LL;
const LL D = 1e4 + 0.5;
LL input() {
    double x;
    scanf("%lf", &x);
    return llround(x * D); // llround()  指的是返回值为 long long, 而参数没有变化
}

功能:将这个数的小数点往右移动 44 位(最后一位要四舍五入),将一个小数变为一个整数,减少失真。(假如想要算出这个数,就将这个数的小数点往左移动 44 位)

那你可能就会问了:假如我就是要算出来根号怎么办?这不难,二分啊!(这里以已知直角边算出斜边为例)

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
LL A, B, C;
int main() {
    scanf("%lld%lld", &A, &B);
    LL l = 0, r = A * A + B * B;
    while (l <= r) {
        LL mid = (l + r) >> 1;
        if (mid * mid <= A * A + B * B) l = mid + 1;
        else r = mid - 1;
    }
    printf("%lld\n", r); // r 是向下取整, l 基本上是向上取整 (l = r + 1) , 但在 r * r = A * A + B * B 的时候,l 就不是向上取整了.
    return 0;
}

这里以一道题为例 ABC191D :给你一个圆的圆心坐标 (X,Y)(X, Y) ,以及圆的半径 RR ,问你园内整点的个数。( X,Y105\vert X \vert ,\vert Y \vert \le 10^5 ,并且 X,YX, Y 小数点后最多有 44 位)

如果用浮点数就会被卡精度,那么我们就这么写:

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
LL read() {
    double x;
    scanf("%lf", &x);
    return llround(x * 10000);
}
LL sq(LL x) { return x * x; }
int main() {
    X = read(); Y = read(); R = read();
    LL ans = 0, RR = sq(R);
    LL high = Y+1, low = Y-1, pos = X - R;
    while (pos % 10000) pos++;
    while (high % 10000) high++;
    while (low % 10000) low--;
    for (; pos <= X + R; pos += 10000) {
        LL A = sq(pos - X);
        while (high >= Y && A + sq(high - Y) > RR) high -= 10000;
        while (high <= Y || A + sq(high - Y) <= RR) high += 10000;
        while (low <= Y && A + sq(low - Y) > RR) low += 10000;
        while (low >= Y || A + sq(low - Y) <= RR) low -= 10000;
        LL add = (high - low) / 10000 - 1;
        if (add > 0) {
            ans += add;
        }
    }
    printf("%lld\n", ans);
    return 0;
}